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ABSTRACT 

The near future of astrophysics involves many large solid-angle, multi-epoch, 

multi-band imaging surveys. These surveys will, at their faint limits, have data 

on large numbers of sources that are too faint to be detected at any individual 

epoch. Here we show that it is possible to measure in multi-epoch data not only 

the fluxes and positions, but also the parallaxes and proper motions of sources 

that are too faint to be detected at any individual epoch. The method involves 

& fitting a model of a moving point source simultaneously to all imaging, taking 

account of the noise and point-spread function in each image. By this method it is 

possible to measure the proper motion of a point source with an uncertainty close 

to the minimum possible uncertainty given the information in the data, which 

is limited by the point-spread function, the distribution of observation times 

(S| (epochs), and the total signal-to-noise in the combined data. We demonstrate 

> 

SDSS Southern Stripe. We show that with our new technique we can use proper 
motions to distinguish very red brown dwarfs from very high-redshift quasars 

qq in these SDSS data, for objects that are inaccessible to traditional techniques, 

and with better fidelity than by multi-band imaging alone. We re-discover all 10 

O known brown dwarfs in our sample and present 9 new candidate brown dwarfs, 

j> identified on the basis of significant proper motion. 

X 

Subject headings: astrometry — methods: statistical — quasars: general — 
stars: kinematics — stars: low-mass, brown dwarfs — techniques: image process- 
ing 



our technique on multi-epoch Sloan Digital Sky Survey (SDSS) imaging of the 
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1. Introduction 

There are many multi-epoch imaging surveys in progress or coming up, which will, 
among other things, deepen our image of the sky and provide information on source variabil- 



ity and proper motions. These surveys include the SDSS Southern Stripe (Abazajian etal 



2009), the Dark Energy Survey, PanSTARRS, LSST, and SNAP. These surveys promise 
proper- motion precisions for well-detected sources on the order of masyr -1 over large parts 
of the sky. For context, a typical halo star at a distance of 10 kpc moving at a transverse 
heliocentric speed of 100 kins" 1 has a proper motion of 2 masyr -1 , and a typical disk star 
at 100 pc and 10 kins" 1 has a proper motion of 20 masyr" 1 . These surveys therefore have 
the capability of revolutionizing our view of the Galaxy and of the Solar neighborhood. 

In most conceptions of a proper-motion measurement, one imagines measuring the po- 
sition of a source in each of several images, taken at different times. A linear trajectory is 
fitted to the positions, relative to some reference frame or set of fixed sources or sources with 
well measured proper motions. In its most straightforward form, this method only works for 
sources bright enough to be detected independently at every epoch — or at least most epochs. 



In a multi-epoch survey like the SDSS Southern Stripe, which has ~ 70 epochs (Abazajian 



etal. 2009), this limits the sources with measured proper motions to a small subset of all 
sources detectable in the combined data, since the combined data reach ~ 2.3 mag fainter 
than any individual epoch; for typical source populations this represents increases in pop- 
ulation size by factors of 5 to 25 at any given signal-to-noise threshold. In this paper we 
present a methodology for measuring in multi-epoch imaging the proper motions of sources 
too faint to detect at any individual epoch. 

There are several different technical regimes for these faint-source proper-motion mea- 
surements. In the "easy" regime, the sources of interest move a distance smaller than or 
comparable to the point-spread function width over the duration of the multi-epoch survey. 
In this regime, the sources are easy to detect in the co-added image, even without taking 
account of their proper motions; proper motions can be determined from processing the 
individual epoch images after detection in the co-added image. There is a "difficult" regime 
in which the sources of interest move substantially more than the width of the point-spread 
function over the duration of the survey. In this regime, the source will not appear at high 
significance in the co-added image if it does not appear at high significance at any epoch, 
because its different appearances in the different individual-epoch images do not overlap. 
In principle, the difficult regime can be addressed by brute force with large computing re- 
sources. In the context of outer Solar-System bodies, brute-force search in the narrow range 



of expected motions is feasible (for example, Bernstein etal. 2004; Fuentes etal. 2008). In 
this paper, we consider only the easy regime. 
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Modeling the data: The traditional method for measuring a stellar proper motion with 
a set of images taken at different times is as follows: Detect the star at each observed epoch; 
measure its centroid (by, for example, finding the peak or first moment of the flux) at each 
observed epoch; and fit a linear motion to the measured positions and times. This procedure 
obtains a proper motion, but it puts an unnecessary requirement on the data: that the star 
be detectable at every epoch. It also puts an unnecessary burden on the data analyst: it 
requires decision making about detection and centroiding of the stars at each epoch, decisions 
that matter at low signal-to-noise, or when faced with data issues such as bad pixels or strong 
variations in noise from pixel to pixel. 

Our new approach is to model all individual-epoch images simultaneously with a single 
point source that is permitted to have a non-zero parallax and proper motion. This approach 
combines the individual-image positional measurement and the determination of the parallax 
and proper motion, and determines all of these simultaneously by making a statistically 
"good" model of the union of all the data. 

In any well-understood imaging survey, each image will have a per-pixel noise model, 
photometric calibration parameters, and a model of the point-spread function. In any suffi- 
ciently small patch of the sky, if the foreground-subtracted intensity in that patch is dom- 
inated by a small number of point sources, it is possible to make an accurate model of all 
of the pixels in the data set that contribute signal to that small patch. In this model of the 
patch, the fluxes, angular positions, parallaxes, and proper motions of the stars in the patch 
are simply parameter values in the well-fitted models. In other words, we are assuming that 
it is possible to model the set of pixels (from all of the images) that contribute to the patch 
with a 6iV-dimensional model that consists of a set of N moving point sources. 

The proper motions determined by image modeling have several advantages over those 
determined by the traditional method: They require fewer decisions about measurement 
techniques (although they do require a good model of the data, including point-spread func- 
tion); they use all of the information in all of the pixels, not just those pixels involved in 
traditional centroiding; they gracefully handle missing data due to bad pixels or cosmic rays 
(assuming the bad pixels have been flagged); they require the investigator to make explicit 
the assumptions about the physical properties of the image and the noise; they can be made 
to properly propagate pixel-value uncertainties into parameter uncertainties (in this case, 
proper motion uncertainties); they are the result of optimization of a well-justified scalar 
objective function (in this case the likelihood). Most importantly for what follows, they 
can be determined in data sets in which the stars are not well detected at any individual 
epoch, but only appear in the combination of the images. In a data set with ~ 70 similar 
epochs (such as the SDSS Southern Stripe), this corresponds to an increase in the number 
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of available targets by factors of 5 to 25 (assuming source populations double to quadruple 
with each magnitude of depth). 

Here we propose, build, test, and use an image-modeling system for the determination 
of stellar proper motions. We show that it can work down to low signal-to-noise ratios and 
that it makes measurements in real data that fully exploit the information available. We also 
use it to discover interesting new astrophysical sources. An approximation to the technique 
used here has been used previously in the Solar System literature (Bernstein etal. 2004). 



Proper-motion and parallax uncertainties: Consider a well-sampled image i with a 
point-spread function of full width at half maximum #fwhm«- The signal-to- noise at which 
the flux of a point source can be measured, [s/n]^ is the sum in quadrature of the signal-to- 
noise contributions from pixels within the point-spread function. A point source measured 
with signal-to-noise [s/n] j in a single image can be centroided with (RMS) uncertainty agj 
of 

<T0,i ~ r / j ; (1) 



details such as the shape of the point-spread function introduce factors of order unity (King 



1983). 



If we have N such images spanning some time interval, we might hope to obtain a proper 
motion estimate with uncertainty limited by the point-spread function, the time interval, 
and the total signal-to-noise 

E [*M 2 (2) 

images i 

in the combination of all the images (we have assumed here that the images i are all in- 
dependent). The relevant time "interval" is not the total time spanned by the data but 
rather S t = ^/Var(i), the standard deviation (root variance) of the times; the best possible 
proper-motion estimates will have uncertainties 

#FWHM /q\ 

where properly 6*fwhm is the square-signal-to-noise weighted mean point-spread function 
full width at half maximum, and S t is the square root of the square-signal-to-noise-weighted 
variance of the times at which the individual epoch images were taken. 

By a similar argument, we hypothesize that the best possible parallax estimates will 
have uncertainties 

[s/n] total 
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where 5\ is the square root of the square-signal-to-noise-weighted variance of the trigono- 
metric functions of the ecliptic longitude A of the Sun (time of year in angle units): 



r2 —— 2 | 2 
°X — °cosA ~r °sinA 



(5) 



Essentially, 8\ describes how well the parallactic ellipse is sampled; an ideal survey for 
parallax measurements will have 5\ ~ 1. 

Disk stars move with respect to one another at velocities of 



30kms (Dehnen & 



Binney 1998 Hogg etal. 2005), that is, on the same order as the velocity of the Earth 



around the Sun. In a multi-epoch survey spanning a small number of years (such as the 
SDSS Southern Stripe), St is of order unity, so for disk stars the parallax and proper motion 
signal-to-noise ratios ought to be comparable in magnitude. However, most surveys sample 
ecliptic longitude A poorly, because of season and scheduling constraints; therefore 5\ is 
usually substantially less than unity, so the signal-to-noise of parallax is smaller than that 
of proper motion. 



2. Method 

The goal is to measure the proper motions and parallaxes of sources detected in multi- 
epoch data. We start with a catalog of detections from a co-addition of the multi-epoch data 
(co-added at zero lag or under an assumption that the sources are static). These detections 
serve as "first-guess" positions for sources in the imaging. We measure the properties of 
these sources by building models of all the individual images, at the pixel level, so that each 
model "predicts" every pixel value in every image at every epoch. 

Some of the candidate sources will not be point sources but rather resolved galaxies, and 
others will not be astronomical sources but will be caused by artificial satellites or imaging 
artifacts. We fit three qualitatively different models, described below. One is of a moving 
point source, one is of an extended galaxy, and one is of a general transient or artifact. 
For each model, "fitting" constitutes optimizing a scalar objective, which is the logarithm 
of the likelihood under the assumption that the per-pixel noise is Gaussian with a known 
variance in each pixel. Under the Gaussian assumption, we can use the different values of 
the log likelihood to perform a hypothesis test based on likelihood ratios. This hypothesis 
test distinguishes point sources from extended galaxies and transients and artifacts. The 
parameters of the best-fitting model are the "measurements" of the source. 

Nothing in what follows fundamentally depends on the assumption of Gaussian noise. 
Data with Poisson errors, for example, can be analyzed the same way but with the objective 
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function changed to the logarithm of the Poisson likelihood. Indeed, any noise model can be 
accomodated, though possibly at the expense of computational simplicity. 

In detail, for each source, we have N small images (patches of what is presumed to be 
a much larger imaging data set) i taken at times U, and we assume that each image has 
reasonable photometric calibration, a noise estimate in each pixel (assumed Gaussian, but 
that could be relaxed in what follows), and correct astrometric calibration or world coordinate 
system (WCS) fixed to an astrometric (RA, Dec) reference frame. From a co-added image 
made from all N single-epoch images we have been given a candidate ( "first-guess" ) position 
(RAj, Decj) for each source j. 

point-source model The first of the three models is that of a point source, moving in 
space and a finite distance from the Solar System. This point source is assumed to have a 
constant flux Sj, a position (RAj, Dec,) at some standard epoch, a parallax ttj and a proper 
motion jlj = (fi a j,fiSj)- In this model and the models to follow, we assume that the sky 
level has been correctly fitted and subtracted from the images, or else that sky errors are 
not strongly covariant with errors in the model parameters. In fitting this model, we find 
the six-dimensional quantity (Sj, RAj, Deo,, 7ij, fj, a j, fisj) that optimizes the scalar objective. 

Given the times U and WCS of the images, any point-source parameter set 
(Sj, RAj, Deo,, TTj, fi a j, H8j), specifies the pixel position of point source j in each image i. 
This position and the (possibly position-dependent) point-spread function model for image 
i permits construction of a pixel-for-pixel model of source j as it ought to appear in image i. 

If we had multi-band imaging (the tests below are on are on single-band images), the 
flux Sj becomes a set of fluxes Skj, one for each bandpass k. In principle, precise fitting 
is complicated by the existence of differential refraction for sources with extreme colors, so 
there are relationships among the fluxes S^, positional offsets, and the airmass or altitude 
of the observations. In the tests below, we are working far enough to the red that there are 
no differential refraction issues at the relevant level of precision. 

Although we have assumed non-varying flux in our model, we should still be able to 
detect and measure moving sources with varying flux. We have not investigated this question, 
but we expect that our method would produce a flux estimate of approximately the mean 
flux measured at the available epochs, and that the point-source model would be preferred 
over the transient model, since the objective function is convex. 

galaxy model Our model of a resolved galaxy is a Gaussian distribution of flux with an 
elliptical covariance parameterized by its radius r,-, eccentricity tj, angle 9j and total flux Sj. 
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For each image, this Gaussian model is convolved with that image's particular point-spread 
function to make a seeing-convolved galaxy model. This seeing-convolved Gaussian galaxy 
model is not a realistic galaxy model, but it is good enough for distinguishing resolved and 
unresolved sources at the faint limit, which is sufficient here. Again, if we had multi-band 
imaging, the flux Sj would be replaced by a set of fluxes Sjk- 

junk model Our model of a transient or imaging artifact is that there is nothing but noise 
in all but one of the images, and that one "junk" image contains many bright pixels. We 
compute this model trivially by computing the chi-squared (x 2 ) contribution for each image 
under the assumption that there is no flux in the image at all. The image with the largest 
X 2 contribution is judged to be the "junk" image and is discarded. In order to keep the 
number of x 2 contributions constant, we replace the "junk" image x 2 by the median of the 
X 2 contributions of the remaining images. 



scalar objective optimization The choices of model, scalar objective, and optimization 
methodology can all be made independently. For the objective function the natural choice 
is the x 2 difference between the model and the data taken over all the pixels that are close 
to the first-guess position in all N images. This objective is analogous to a logarithm of 
a likelihood ratio; it is exact if the noise in the image pixels is Gaussian and independent, 
with known variances (which can vary from pixel to pixel). For optimizing this objective 
function, we use the Levenberg-Marquardt method (Levenberg 1944 Marquardt 1963). 



hypothesis test In the approximation that the noise is Gaussian, the best fits for each 
of the three models can be compared via the best-fit values of the x 2 scalar objective. If 
the three models are equally likely a priori and if they have the same number of degrees of 
freedom, then one model is confidently preferred over another if it has a best-fit x 2 value 
smaller by an amount Ax 2 ^> 1. Of course the models are not equally likely a priori, but for 
for the vast majority of sources, the differences in x 2 are so large that no reasonable prior 
would change the results of our hypothesis test. 

Note that there is some degeneracy in our models: a galaxy model with zero radius and 
a star model with zero proper motion and parallax produce exactly the same predictions, and 
thus our hypothesis test cannot distinguish them. This could be remedied by placing prior 
probabilies over the model parameters — for example, penalizing tiny galaxies — but since we 
are not concerned with the region of parameter space where this occurs, we have not done 
this. 
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Rather than explicitly including a junk model, we could instead place a threshold on 
the likelihood of the star and galaxy models: junk data will be poorly fit by the star and 
galaxy models and thus will have tiny likelihood. In general we have found that image sets 
for which the junk model is preferred clearly contain artifacts or transients; the method is 
not sensitive to the details of the junk model. 

jackknife error analysis In principle, the region in parameter space around the best-fit 
point where \ 2 is within unity of the minimum provides an estimate of the uncertainties 
in the fitted parameters. However, this estimate is only good when the model is a good 
fit; many error contributions in real data come from source variability, poorly known data 
properties (such as pixel uncertainty or point-spread-function estimates that are in error) 
and unflagged artifacts in the data. For this reason, we use (and advocate) a "jackknife" 
technique for error analysis. 

The jackknife technique is to perform the analysis on the N subsets of the N images 
created by leaving one image out. The complete fit of the three models is performed on each 
of the N leave-one-out subsets and parameters are measured. The uncertainty estimate a p 
for any fitted parameter p is related to the N leave-one-out measurements pi (made leaving 
out image i) by 



where (p) is the mean of the leave-one-out measurements pi. The jackknife technique auto- 
matically marginalizes the error estimates over the other parameters, and provides a prop- 
erly marginalized estimate of any multi-parameter covariance matrix by the generalization 
of equation ^ in which the square is changed into the d x d matrix outer product of the 
"vectors" made from the d parameters for which the covariance matrix is desired. Of course 
when d is large, the jackknife will not accurately sample all degrees of freedom available in the 
covariance matrix, but provided N is large enough, it will sample the dominant eigenvectors 
(the principal components). 

Implementation notes Our code is implemented in Python and uses the Django web 
framework, which provides powerful database and web server integration. This allows us to 
quickly and easily manage and visualize the data and results. Combined with the scientific 
data analysis packages scipy and numpy and the plotting package matplotlib, this yielded 
a powerful software development environment. 
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For optimization, we use the Levenberg-Marquardt implementation levmar (version 2.2; 
Lourakis 2004) with Python bindings pylevmar (revision 313; Tse 2008). In this Python 
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environment, analysis takes on the order of seconds for each source (30 epochs, 15 x 15 
images), but this could be sped up substantially by implementing some of the core operations 
in C. 



Tests on real data 



For test data, we make use of the SDSS Southern Stripe (SDSSSS), a multi-epoch survey 



undertaken as part of SDSS-II (Adelman-McCarthy etal. 2008 Abazajian etal. 2009). The 



SDSSSS data are part of The Sloan Digital Sky Survey (Gunn etal. 1998 York etal. 2000); 
it involves ugriz CCD imaging of ~ 250 deg 2 on the Equator in the southern Galactic cap. 



All the SDSSSS data processing, including astrometry (Pier etal. 2003), source identifica- 



tion, deblending and photometry (Lupton etal. 2001), and calibration (Smith etal. 2002 



Padmanabhan etal. 2008) are performed with automated SDSS software. 



The SDSSSS data have been found to have a small astrometric drift (Bramich etal 



2008), because astrometric calibration was performed at a single, slightly inappropriate epoch 
(Pier etal. 2003). This drift, for which we are making no correction, is at the 10 masyr -1 



level; at the precision of this study it does not change any of the conclusions below. 

In general, the hypothesis test we perform requires that the variance of the noise be 
properly estimated on a pixel-by-pixel basis. These are based on an SDSS imaging noise 
model, with the adjustment that pixels that have been corrupted by cosmic rays or other 
defects are given infinite variances (vanishing contribution to x 2 )- Occasionally there are 
unidentified cosmic rays in the data. These lead to localized regions with very large contri- 
butions to x 2 - When one of these noise defects appears in the data near one of the targets, 
it sometimes causes a source which is truly a galaxy or a star to be assigned "junk" status. 
After by-eye inspection of cutouts, we estimate this rate to be on the order of < 0.5 percent 
for this data source; the rate of such problems increases with the number of epochs and the 
image cut-out size (the total number of pixels in the fit). 



For some of the sources we have UKIRT Infrared Deep Sky Survey (UKIDSS; Lawrence 



etal. 2007) data. UKIDSS uses the UKIRT Wide Field Camera (Casali etal. 2007) with 



the infrared photometric system described by Hewett etal. (2006), and automated data 



processing and archiving (Irwin etal. 2008 Hambly etal. 2008). The UKIDSS data used 



here comes from the fourth data release. 

Very red point sources in deep optical imaging — for example, ^-band-only sources in the 
multi-epoch SDSS Southern Stripe — include both very cool dwarfs and very high redshift 
quasars. In principle these can be distinguished with parallax and proper-motion estimates. 
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For this reason, we performed a test on z-only point sources in the SDSS Southern Stripe. 
The parent sample is point sources from the SDSSSS "Co-add Catalog" (J. Annis etai, in 
preparation) that have [i — z] > 2 mag and [r — z] > 2 mag. This criterion selects quasars at 



5-8 < z < 6.5 as well as cool dwarfs with spectral types ranging from mid-L to T (Fan etal. 



2001 and references therein). Hotter brown dwarfs, stars, and lower-redshift quasars have 
significant emission in the i band, giving them bluer i — z colours, while the emission features 
of cooler dwarfs and higher-redshift quasars lie mostly redwards of the SDSS z bandpass. 



The Co-add Catalog uses the asinh magnitude (or "luptitude") scale (Lupton etal 



1999), so it is possible to select objects based on color even for objects that are not detected 
in one of the bands. The version of the Co-add Catalog we are using is from SDSS Data 
Release 7 (DR7), and includes 20 to 40 epochs over ~ 250 deg 2 . 

Since we are interested in distinguishing cool brown dwarfs from high-redshift quasars, 
we require <r M to be small (equation [3J. We therefore cut our parent sample to have 
z < 21 mag, leaving roughly 150 sources. This cut allows us to reach, with moderate 
signal-to- noise, slightly fainter sources than are detectable in the single-epoch images. In a 
future paper, we plan to relax this cut, which should yield considerably more brown dwarf 
candidates at smaller signal-to-noise levels. Of our 150 parent candidates, some turn out to 
be caused by an imaging artifact or transient in one of the N epochs, and some turn out 
to be galaxies or stars with mis-measured colors because of data artifacts or inaccuracies in 
deb lending nearby objects. 

Each of the catalog sources has a nominal position and a z-band magnitude in the Co- 
add Catalog. For each 2-only source, we cut out 15 x 15 pix 2 patches of every SDSS image 
at the nominal position. For each tiny image, we construct a tiny local world-coordinate- 
system description of the astrometric calibration of that patch using the SDSS pipeline 
astrometric calibration. We subtract the local value of a smoothly fit sky level (M. Blan- 
ton, in preparation) and multiply each tiny image by a constant, based on the pipeline 
calibration information, to place it on a common photometric calibration scale in intensity 
units (energy per unit solid angle per unit area per unit time per unit frequency). To the 
SDSS pipeline-reconstructed point-spread function (PSF) in each tiny image we fit a single- 
Gaussian approximate model, which is not a good fit to the PSF at high precision, but which 
is sufficient for modeling sources at low signal-to-noise. 

In this work, we create the 15 x 15 pix 2 cutouts from the same 20 to 40 epochs that 
are used in the Co-add Catalog; in future work we plan to use the ~ 70 epochs that have 
become available in DR7. 

We chose 15 x 15 pix 2 patches so that a source with proper motion of ~ 0.5 arcsecyr -1 
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would remain in the patch. Our method still works if the source leaves the patch — indeed, 
our fastest-moving candidate does this — but we gain no information from the epochs in which 
the source has left the patch, so the signal-to-noise of our parameter estimates will be less 
than optimal when this happens. We could choose to use larger cutouts; the only difficulty is 
that if the patch contains more than one source, our model will try to explain the brightest 
source (because this will decrease x 2 the most). This could perhaps be remedied by adding 
a prior on the source position, but since the SDSSSS data are from well below the galactic 
equator, stellar density is low and we have not found this to be necessary. Alternatively, in 
cases where the source leaves the original patch we could produce new cutouts that track 
the source motion. 

Figures [I] through [4] illustrate our approach by showing the results of the (moving) point- 
source and (static) galaxy model fits to four sources in the the SDSSSS data. In these figures 
we show all the individual 15 x 15 images from the individual epochs, and the best-fit point- 
source and galaxy parameters. In these figures, we visualize the distribution of acceptable 
parameters around the best-fit values through sampling. We also show mean images and 
mean residual maps in the static and moving coordinate systems. These figures demonstrate 
heuristically that the hypothesis test is effective at separating sources of different types, even 
when the source is not apparent at high signal-to-noise at any individual epoch. 



In Figure [5] we show the overall results from application of our techniques to the [i — z] > 
2 mag sources in the SDSSSS: We show proper-motion measurements and jackknife estimates 
of our uncertainties as a function of z-band magnitude. Known quasars and brown dwarfs 
are marked. Our measurements clearly separate the known quasars and brown dwarfs on 
the basis of proper motion alone. All known brown dwarfs in the sample obtain significant 
non-zero proper motion measurements, and all known high-redshift quasars in our sample 
obtain proper motion measurements consistent with zero. The sources in our sample that 
have significant motions and have not been previously identified as brown dwarfs are our 
new brown dwarf candidates. In Figure [6] we show the UKIDSS and SDSS [z — J] colors of 
the sources for which we have UKIDSS J measurements, with the known brown dwarfs and 
quasars and our new brown dwarf candidates marked. Many of the sources in these Figures 
are undetectable (or not detectable reliably) at individual epochs; the single-epoch 5-sigma 



detection limit is roughly z = 20.5 mag in good seeing conditions (Abazajian etal. 2009). 



In Figure [7| our jackknife estimates of our measurement uncertainties are compared 
to approximate estimates of the total information content in each source's data set, made 
with an approximation to equation (J3J). If our uncertainty estimates are correct (as we 
demonstrate that they are, below), this shows that we come close to attaining the accuracy 
available. 
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4. Tests on artificial data 

To demonstrate that our jackknife error estimates are reasonable, and that our code is 
optimizing the models correctly, we performed some tests on synthetic data. We selected a 
subset of the SDSSSS candidate objects for which we found reasonable fits to a moving point 
source model. For each candidate, we generated a stack of images by generating, for each 
image in the original stack, the image predicted by our point-source model, given the WCS, 
point-spread function, time, and noise amplitude of the image. This is a good test set because 
it has the same imaging properties as the original data and the same distribution of point- 
source parameters as the sources we want to be able to discover. Since the synthetic images 
are generated using our image model, this test shows how our algorithm would perform if 
our modeling assumptions were exactly correct. 

After running our optimization code on these synthetic images, we compare our errors — 
the differences between the true and estimated moving-point-source parameters — to the jack- 
knife estimates of our uncertainties. In Figure [8] we show that the errors are consistent with 
the uncertainty estimates. This shows that when our assumptions about the data are correct, 
we do measure the proper motions as accurately as our jackknife errors indicate. 



5. Discussion 

We have shown that straightforward image modeling permits the measurement of ap- 
parent motions, especially the proper motion and parallax of a source in multi-epoch data, 
even when the source is too faint to be reliably detected or centroided at any individual 
epoch. The results of this project are not surprising; indeed what is surprising is how rarely 
the measurements of stellar motions are made by comprehensive data modeling. 

We demonstrated the technique on real and artificial data. In the process of perform- 
ing these tests we showed that spectrosopically confirmed quasars and brown dwarfs can be 
perfectly distinguished with proper motions measured by this technique. Working without 
proper motions, but with Co-add Catalog sources and a significant amount of near-infrared 
imaging follow-up, a group has followed up the z-only sources most likely to be high-redshift 



Chiu etal. 


2008 




Jiang etal. 


2008) 



This project, even after infrared imaging, 



found — after expensive spectroscopic follow-up — that some of the high-redshift quasar can- 
didates selected on the basis of visible and near-infrared imaging are in fact nearby brown 
dwarfs. We have shown that all of these spectroscopically confirmed brown dwarfs have sig- 
nificantly measured (> 5 sigma) non-zero proper motions by the technique shown here (and 
are reported in Table [11. None of the spectroscopically confirmed high-redshift quasars do. 
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Use of this technique could have been used to substantially increase the efficiency of either 
quasar or brown-dwarf searches in this data set. 

In performing this demonstration, we have independently identified all 10 known brown 



dwarfs (Fan etal. 2000 Geballe etal. 2002; Hawley etal. 2002 Berriman etal. 2003 Knapp 



eta/.|2004||Chiu eta/.|2008||Metchev etal.\2008\) in our parent sample, and we have discovered 
9 new candidate brown dwarfs, presented in Table [TJ Based on our analysis, these objects 
have a high probability of being brown dwarfs. It would be desirable to separate disk 
dwarfs from halo dwarfs — the fastest angular movers tend to be halo members (for example, 
Lepine etal. 2003) — but the time cadence of the SDSSSS data is such that parallaxes are 
not measured well. Two of the dwarfs we rediscover — 2MASS J010752. 42+004156. 3 and 



2MASS J020742.84+000056.4— have previously measured parallaxes flVrba eta/.||2004[ ); the 
measurements are consistent with our upper limits. 

Our tests show that the uncertainty in the proper-motion measurement made by image 
modeling is consistent with the best possible uncertainties given the angular resolution and 
photometric sensitivity of the combination of all images in the multi-epoch data set. These 
tests effectively show that such measurements can be made for objects that are fainter than 
those available to traditional methods that require source detection at every epoch. In 
imaging with N equally sensitive epochs, we are able to measure objects that are fainter by 
Am magnitudes: 



Am = -bg 2J ([s/n] i )+log 2 .5([s/w]totai) 
= log 2 . 5 (ViV) 
~ 0.55 logiV mag . 



(7) 



(9) 



This advantage amounts to 1 mag for surveys with 6 similar epochs, and 1.6 to 2.0 mag in 
data with 20 to 40 epochs (such as the data used here). In the ~ 70 epochs available in SDSS 
DR7, it reaches 2.3 mag. Several of the high-redshift quasars and brown dwarfs analyzed in 
this study were only detectable in the combination of all of the multi-epoch images. 

The depth advantage of image modeling is most dramatic in surveys with very large 
numbers of epochs, as is expected for LSST. In general the number of interesting sources is 
a strong function of depth (factors of 2 to 4 per magnitude), so the "reach" of the image- 
modeling technique is a strong function of the number of epochs. 

One limitation of the work presented here is that we used the zero-proper-motion image 
"stack" for source detection and therefore will only have in the candidate list objects with 
small proper motions. Faint stars and dwarfs with proper motions large enough that they 
move the width of the PSF between epochs, or some significant fraction of that, are harder 
to find, because they don't appear in the stack at much higher signal-to-noise than they 
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appear in any individual-epoch image. In future work we hope to address the detection 
and measurement of these fast-moving but very faint sources. Approximations have been 
executed in the search for Solar System bodies (for example, Bernstein et a/.|2004 ). Certainly 
a reliable system for discovery in this regime would have a big impact on future surveys like 
PanSTARRS and LSST. 
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Name 


R A 


Dec 


flux 


parallax 


dRA/dt 


dDec/dt 


notes 




dog 


dcg 


s/n 


arcsec 


arcsec yr — 1 


arcsec yr - 1 








SDSS J001608. 47-004302. 9 


4.03537 


-0.71733 


46.4 


+ 0.007 ± 0.058 


+0.134 ± 0.015 


-0.020 ± 0.007 


2 


u 


BD 


SDSS J005212. 29+001216.0 


13.05131 


+0.20465 


29.6 


-0.053 ± 0.069 


-0.165 ± 0.014 


-0.211 ± 0.010 


2 


u 


BD 


SDSS J010407.68-005329.1 


16.03195 


-0.89128 


59.1 


-0.041 ± 0.068 


+0.460 ± 0.010 


-0.017 ± 0.008 


2 


u 


BD 


SDSS .1010752.59+004156.0 


16.96899 


+0.69905 


59.0 


+0.044 ± 0.070 


+0.644 ± 0.014 


+0.085 ± 0.013 


2 




BD 


SDSS J020333. 28-010813.1 


30.88881 


-1.13683 


18.6 


-0.019 ± 0.063 


+0.354 ± 0.019 


-0.005 ± 0.012 




u 


BD 


SDSS J020742. 85+000055. 6 


31.92867 


+0.01561 


27.8 


+0.078 ± 0.072 


+0.163 ± 0.017 


-0.029 ± 0.010 


2 


u 


BD 


SDSS J023617. 95+004853. 5 


39.07492 


+0.81501 


74.4 


+0.025 ± 0.021 


+0.134 ± 0.006 


-0.166 ± 0.005 


2 


u 


BD 


SDSS J033035. 23-002537. 2 


52.64687 


-0.42678 


76.5 


+0.040 ± 0.026 


+0.390 ± 0.007 


-0.360 ± 0.008 


2 


u 


BD 


SDSS .1214046.48+011258.2 


325.19385 


+ 1.21633 


28.3 


-0.064 ± 0.143 


-0.085 ± 0.012 


-0.215 ± 0.007 


2 




BD 


SDSS J224953. 45+004403. 9 


342.47285 


+0.73458 


44.4 


-0.003 ± 0.142 


+0.084 ± 0.010 


+0.011 ± 0.009 


2 


u 


BD 


SDSS J001836. 46-002559. 9 


4.65204 


-0.43315 


27.6 


+0.038 ± 0.137 


+0.179 ± 0.019 


-0.029 ± 0.017 




I' 




SDSS J011014. 40+010618. 5 


17.55990 


+1.10534 


32.0 


-0.017 ± 0.089 


+0.542 ± 0.023 


+0.013 ± 0.015 




u 




SDSS J011417.92-003437.9 


18.57481 


-0.57704 


19.4 


+0.090 ± 0.064 


-0.093 ± 0.021 


-0.077 ± 0.018 




u 




SDSS J021642. 94+004005.1 


34.17907 


+0.66828 


52.0 


-0.021 ± 0.023 


-0.069 ± 0.011 


-0.093 ± 0.009 




I' 




SDSS J023047.97-002600.4 


37.69996 


-0.43332 


29.9 


-0.012 ± 0.045 


+0.127 ± 0.008 


-0.003 ± 0.010 




I' 




SDSS J215919. 95+003309.0 


329.83326 


+0.55260 


21.5 


+0.120 ± 0.158 


+0.155 ± 0.025 


+0.100 ± 0.018 








SDSS J234730. 64-002912.0 


356.87782 


-0.48653 


14.7 


+0.095 ± 0.120 


-0.082 ± 0.020 


-0.090 ± 0.026 




u 




SDSS J234841. 38-004022. 9 


357.17250 


-0.67289 


39.7 


+0.161 ± 0.084 


+0.097 ± 0.025 


-0.125 ± 0.035 




I' 




SDSS J235410. 42+004315. 9 


358.54362 


+0.72131 


59.6 


-0.061 ± 0.077 


+0.053 ± 0.013 


-0.063 ± 0.010 




I' 




Table 1: Well-fit [i 


- z] > 2 


mag sources 


in the SDSS 


Southern Stripe with 


proper 



> 60 masyr -1 measured at high confidence (> 3.5 sigma). RA, Dec positions have equinox 
J2000.0 but are computed for MJD 53000. The note "2" indicates that there is a nearby 
entry in the 2MASS point-source catalog ( Skrutskie eta/.|2006 ), while "U" indicates a nearby 



source in the UKIDSS catalog (Lawrence etal. 2007). "BD" indicates objects that are 



spectroscopically-confirmed brown dwarfs (Fan etal. 2000 Geballe etal. 2002 Hawley etal. 



2002 Berriman etal. 2003 Knapp etal. 2004 Chiu etal. 2008; Metchev etal. 2008). The 



9 sources in the lower part of the table (not marked with "BD' : 
candidates. 



are new brown-dwarf 
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Name 


SDSS u 


SDSS g 




SDSS r 


SDSS i 


SDSS z 






mag 


mag 




mag 


mag 


mag 


SDSS 


J001608. 47-004302. 9 


24.26 ± 0.29 


25.13 ± 


16 


23.88 ± 0.18 


21.31 ± 0.03 


19.23 ± 0.01 


SDSS 


J005212. 29+001216.0 


24.35 ± 0.26 


25.03 ± 


13 


23.88 ± 0.10 


21.92 ± 0.03 


19.62 ± 0.02 


SDSS 


J010407.68-005329.1 


25.02 ± 0.23 


25.17 ± 


22 


24.12 ± 0.13 


22.07 ± 0.04 


19.83 ± 0.02 


SDSS 


J010752. 59+004156.0 


24.31 ± 0.21 


25.17 ± 


13 


24.09 ± 0.11 


21.87 ± 0.03 


19.11 ± 0.01 


SDSS 


J020333. 28-010813.1 


24.60 ± 0.24 


25.18 ± 


21 


24.60 ± 0.14 


23.91 ± 0.14 


20.85 ± 0.04 


SDSS 


J020742. 85+000055. 6 


24.43 ± 0.25 


25.25 ± 


13 


25.04 ± 0.13 


24.15 ± 0.13 


20.43 ± 0.04 


SDSS 


J023617. 95+004853. 5 


24.83 ± 0.20 


24.89 ± 


1L> 


23.85 ± 0.09 


21.69 ± 0.03 


19.01 ± 0.01 


SDSS 


J033035. 23-002537. 2 


24.66 ± 0.25 


25.01 ± 


13 


23.06 ± 0.05 


20.88 ± 0.01 


18.79 ± 0.01 


SDSS 


J214046. 48+011258. 2 


24.60 ± 0.29 


24.67 ± 


15 


23.55 ± 0.11 


21.16 ± 0.02 


19.11 ± 0.01 


SDSS 


J224953. 45+004403. 9 


24.01 ± 0.21 


25.24 ± 


15 


23.79 ± 0.11 


21.61 ± 0.03 


19.53 ± 0.01 



SDSS J001836. 46-002559. 9 24.65 + 0.27 25.09 + 0.14 24.29 + 0.19 23.23 + 0.11 20.43 + 0.03 

SDSS J011014. 40+010618. 5 24.95 + 0.25 25.27 + 0.13 24.28 + 0.13 22.29 + 0.05 20.03 + 0.02 

SDSS J011417. 92 — 003437.9 24.55 + 0.32 25.10 + 0.13 24.76 + 0.15 22.46 + 0.09 20.40 + 0.03 

SDSS J021642. 94+004005.1 24.82 + 0.20 24.98 + 0.12 24.14 + 0.11 22.16 + 0.04 20.01 + 0.02 

SDSS J023047. 97— 002600.4 24.52 + 0.25 24.89 + 0.12 24.13 + 0.11 22.36 + 0.04 20.22 + 0.02 

SDSS J215919. 95+003309.0 24.72 + 0.25 24.89 + 0.15 24.20 + 0.15 22.59 + 0.09 20.58 + 0.04 

SDSS J234730. 64-002912.0 24.46 + 0.28 25.03 + 0.14 24.78 + 0.15 22.93 + 0.08 20.90 + 0.05 

SDSS J234841. 38-004022. 9 24.76 + 0.29 25.07 + 0.16 23.93 + 0.12 21.95 + 0.04 19.93 + 0.02 

SDSS J235410. 42+004315. 9 24.04 + 0.19 25.03 + 0.22 23.24 + 0.10 21.17 + 0.07 19.11 + 0.01 



Name 


UKIDSS y 


2MASS J 


UKIDSS Jl 


2MASS H 


UKIDSS H 


2MASS K s 


UKIDSS K 




mag 


mag 




mag 


mag 


mag 


mag 


mag 


SDSS J001608. 47-004302. 9 


17.69 ± 0.02 


16.33 ± 


12 


16.30 ± 0.01 


15.23 ± 0.11 


15.35 ± 0.01 


14.54 ± 0.09 


14.50 ± 0.01 


SDSS .1005212.29+001216.0 


17.87 ± 0.03 


16.36 ± 


11 


16.54 ± 0.02 


15.56 ± 0.13 


15.81 ± 0.01 


15.46 ± 0.16 


15.21 ± 0.01 


SDSS J010407.68-005329.1 


17.87 ± 0.03 


16.53 ± 


13 


16.54 ± 0.01 


15.64 ± 0.14 


15.80 ± 0.02 


15.33 ± 0.17 


15.24 ± 0.02 


SDSS .1010752.59+004156.0 




15.82 ± 


06 




14.51 ± 0.04 




13.71 ± 0.04 




SDSS J020333. 28-010813.1 


18.99 ± 0.08 






17.70 ± 0.04 




16.89 ± 0.03 




16.26 ± 0.03 


SDSS J020742. 85+000055. 6 


18.03 ± 0.03 


16.80 ± 


16 


16.74 ± 0.01 


16.40 


16.81 ± 0.04 


15.41 


16.73 ± 0.05 


SDSS .1023617.95+004853.5 


17.18 ± 0.02 


16.10 ± 


08 




15.27 ± 0.07 


15.14 ± 0.01 


14.67 ± 0.09 


14.56 ± 0.01 


SDSS J033035. 23-002537. 2 


16.50 ± 0.01 


15.31 ± 


05 


15.21 ± 0.00 


14.42 ± 0.04 


14.49 ± 0.00 


13.84 ± 0.05 


13.77 ± 0.00 


SDSS .1214046.48+011258.2 




15.89 ± 


08 




15.31 ± 0.09 




14.42 ± 0.08 




SDSS .1224953.45+004403.9 


17.78 ± 0.03 


16.59 ± 


12 


16.45 ± 0.01 


15.42 ± 0.11 


15.33 ± 0.01 


14.36 ± 0.07 


14.41 ± 0.01 


SDSS JO01836. 46-002559. 9 


18.73 ± 0.08 






17.67 ± 0.05 




16.61 ± 0.04 






SDSS J011014. 40+010618. 5 








16.90 ± 0.02 




16.34 ± 0.03 




15.73 ± 0.02 


SDSS J011417.92-003437.9 


18.85 ± 0.08 






17.89 ± 0.04 




17.43 ± 0.07 




16.73 ± 0.06 


SDSS J021642. 94+004005.1 


18.41 ± 0.05 






17.31 ± 0.03 




16.97 ± 0.04 




16.51 ± 0.04 


SDSS J023047.97-002600.4 








17.35 ± 0.03 




16.64 ± 0.04 




16.08 ± 0.04 


SDSS J215919. 95+003309.0 


















SDSS J234730. 64-002912.0 


19.56 ± 0.11 






18.37 ± 0.08 




17.87 ± 0.09 




17.18 ± 0.08 


SDSS J234841. 38-004022. 9 








17.50 ± 0.03 




16.89 ± 0.03 




16.40 ± 0.04 


SDSS J235410. 42+004315. 9 


17.82 ± 0.03 






17.01 ± 0.02 




16.38 ± 0.03 




15.96 ± 0.03 



Table 2: Photometric properties for the sources in Table [TJ by association with sources in 
SDSS, 2MASS, and UKIDSS. SDSS magnitudes are given in the SDSS system, which is 
very close to the AB system (See |http : //www, sdss . org/dr7/al gorithms/f luxcal .html# 
sdss2ab); UKIDSS magnitudes are given as retrieved from the WFCAM science archive, 
i.e., they are Vega-based. 




Fig. 1. — [caption on next page] 
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Fig. 1. — [on previous page] The results of fitting the SDSS Southern Stripe multi-epoch 
imaging data on SDSS J020333. 28-010813.1, a spectroscopically confirmed brown dwarf 
(Knapp etal. 2004, L. Jiang, private communication) and [i — z] > 2 mag source. The 
top set of panels — labeled by observation MJD — show the individual epoch 15 x 15-pixel 
sub-images; note that the source is not clearly detectable at every epoch. The middle di- 
agrams show the output of fitting a moving source (left panel) or a resolved galaxy (right 
panel). On the moving-source diagram, the best fit path of the moving point source is 
shown as a thick black line, the thinner grey lines show alternative paths sampled from 
the jackknife-inf erred posterior distribution of trajectories consistent with the data; that is, 
the grey lines effectively show the uncertainty interval. The thick black and thin grey lines 
contain wiggles with a period of one year (or the pixel distance of one year at that path's 
proper motion) because each is the realization of a trajectory with finite proper motion and 
parallax. It can be seen from this panel that this source has a well-measured proper motion 
but not a well-measured parallax, because the grey lines do not share a common parallax. 
The PSF FWHM sizes of the individual epoch images are shown as circles centered on the 
positions the point source would have on the best-fit path. On the galaxy diagram, the 
mean-PSF-convolved galaxy model is shown as a black ellipse, and the grey ellipses sample 
the jackknife- inferred distribution of galaxy models consistent with the data. Note that for 
this source, the point-source model is a much better fit than the galaxy model (with y 2 
difference 169), so the point-source model is favored. Along the bottom, the leftmost panel 
(data stack) shows the data co-added (weighted by per-pixel inverse variance). The second 
and third panels (star stack and galaxy stack) show the star and galaxy models co-added at 
zero lag (no proper motion compensation). The fourth panel (star resid) shows the residuals 
(data minus model) for the point-source model, and the fifth panel (galaxy resid) the resid- 
uals for the galaxy model, both co-added at zero lag. The sixth panel (data shifted) shows 
the data co-added with the best-fit proper motion compensated. The seventh (star shifted) 
shows the point-source model, co-added with the best-fit proper motion compensated. The 
final panel (shifted resid) shows the residuals for the point-source model co-added with the 
best-fit proper motion compensated. The motion-compensating shifts are rounded to the 
nearest pixel. Note that these co-added figures are only shown for purposes of illustration; 
our method never co-adds images. 
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Fig. 2.— The same as Figure [T] but for SDSS J020332.35+001228.6, a spectroscopically 
confirmed z ~ 6 quasar (Jiang etal. 2008) and [i — z] > 2 mag source. Here the point- 



source model is favored, but the inferred proper motion (best-fit value or any sample from 
the distribution) is very small; the wavy paths each span ~ 100 yr in time; they have 
different position angles because when the magnitude of the proper motion is constrained to 
be near zero, the direction is not well constrained. The galaxy model is disfavored by a small 
amount. The amount is small because the best-fit galaxy model is a non- moving compact 
source, which is not dissimilar to the nearly non-moving best-fit point source. 




Fig. 3. — The same as Figure [T] but for a faint galaxy. Here the galaxy model is favored. 
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Fig. 4. — The same as Figure [T] but for a spurious source caused by a blinking artificial 
satellite. Here the "junk" model is favored over the point-source and galaxy models. Note 
that the jackknife errors in the point-source and galaxy models are very large. This source 
actually can be removed by flag-checking in the SDSSSS Coadd Catalog, but we show it here 
for illustrative purposes. 
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Fig. 5. — Proper motion magnitude (angular speed) as a function of SDSSSS Coadd Catalog 
z-band magnitude for [i — z] > 2 mag sources in the SDSS Southern Stripe that are pref- 
erentially described as point sources (by our x 2 hypothesis test). The uncertainty regions 
are shown as transparent ellipses. The spectroscopically confirmed high-redshift quasars 
d Jiang e£al.||2008D and brown dwarfs ( [Fan et aZ.|[2000| |Geballe et al|2002| |Hawley etal\\2QQ2 



Berriman etal. 2003 Knapp etal. 2004; Chiu etal. 2008 Metchev etal. 2008) are shown in 



color. Every one of the brown dwarfs has a significantly measured proper motion; none of 
the quasars do. Other brown-dwarf candidates are clearly visible as significant movers (see 
also Table [TJ. Note that the single-epoch detection limit is approximately z = 20.5 mag in 
good seeing conditions. 
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Fig. 6. — UKIDSS and SDSSSS [z — J] color plotted against SDSSSS z-band magnitude, for 
the sources in Figure [5] that are detected in the UKIDSS J band. All of the significantly- 
moving objects have the very red [z — J] colors of brown dwarfs; the likely brown dwarfs 
could have been identified by their proper motions and SDSSSS Coadd Catalog colors alone. 
Note that the z magnitude is in the AB system, while the J magnitude is Vega-based. 
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Fig. 7. — Comparison of the jackknife-estimated proper-motion uncertainties to "best-case" 
values estimated from general principles of information in the imaging: The information es- 
timate is the mean (square-signal-to-noise-weighted) imaging point-spread function FWHM 
divided by the signal-to-noise of the flux measurement (taken to be a proxy for the total 
detection signal-to-noise), divided by the root- variance of the time span. This figure shows 
that the measurements are roughly as precise as they can be, given the information content 
of the set of images; see equation (|3]). 
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Fig. 8. — Accuracy of jackknife uncertainty estimates for fits to artificial data made with 
known point-source properties. Errors in fit parameters (fit minus true) have been divided 
by jackknife uncertainties. The artificial data sets have identical imaging properties (noise 
amplitude, WCS, and PSF) to SDSS Southern Stripe sources, but contain artifical images 
made with point sources with true positions, fluxes and motions derived from the data as 
described in the text. The cloud of points is centered near (0, 0), is circularly symmetric, and 
appears roughly Gaussian: our estimates are unbiased, uncorrected, and have the expected 
magnitude and distribution of error. 



